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Abstract 

Diffusion MRI is a well established imaging modality providing a powerful 
and innovative way to non-invasively probe the structure of the white matter. 
Despite the potential of the technique, the intrinsic long scan times of these 
sequences have hampered their use in clinical practice. For this reason, a wide 
variety of methods have been proposed recently with the aim to shorten acqui- 
sition times. Among them, spherical deconvolution approaches have gained a 
lot of interest for their ability to reliably recover the intra- voxel fiber configura- 
tion with a relatively small number of data samples. To overcome the intrinsic 
instabilities in solving the deconvolution problem, these methods make use of 
regularization schemes generally based on the assumption that the fiber orienta- 
tion distribution (FOD) to be recovered is sparse, either explicitly or implicitly. 
In particular, convex optimization methods have recently been advocated in 
a compressed sensing perspective for FOD reconstruction from accelerated ac- 
quisitions. In this paper, we propose to exploit further the versatility of this 
powerful framework with the aim to exploit sparsity more optimally. We de- 
fine a new convex minimization problem for FOD reconstruction through a 
constrained formulation between sparsity prior and data, also making use of a 
reweighting scheme. The method has been tested on both synthetic and real 
data. Experimental results indicate that this approach provides more robust 
and accurate estimates than the state of the art in terms of both the number 
and orientation of fiber compartments in each voxel. 
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1. Introduction 

Fiber-tracking is probably one of the most fascinating applications in dif- 
fusion MRI, gathering a lot of attention since the beginning for its ability to 
reconstruct from the acquired data the main neuronal bundles of the brain. In 
fact, the diffusion of the molecules in the white matter can be exploited for 
mapping the brain connectivity, and structures otherwise invisible with other 
imaging modalities can be highlighted. The study of this structural connectivity 
is of major importance in a fundamental neuroscience perspective, for developing 
our understanding of the brain, but also in a clinical perspective, with particular 
applications for the study of a wide range of neurological disorders. 

The most powerful acquisition modality is diffusion spectrum imaging (DSI) 
(Wedeen et al., 2005). It relies on cartesian signal sampling and is known to 
provide good imaging quality, but it is too time-consuming to be of real in- 
terest in a clinical perspective. Diffusion tensor imaging (DTI) (Basser et al., 
1994) is always preferred instead. DTI is a very fast model-based technique 
providing valuable diagnostic information but, on the contrary, it is unable to 
model multiple fiber populations in a voxel. In a global connectivity analysis 
perspective, this constitutes a key limiting factor. Accelerated acquisitions re- 
lying on a smaller number of samples while providing accurate estimations of 
the intra-voxel fiber configuration thus represent an important challenge. 

Recently, an increasing number of high angular resolution diffusion imag- 
ing (HARDI) approaches have been proposed for tackling this problem. In 
particular, spherical deconvolution (SD) based methods formed a very active 
area in this field (Tournicr et al., 2004; Alexander, 2005; Tournier et al., 2007; 
Dcll'acqua et al., 2007). These methods rely on the assumption that the signal 
attenuation acquired with diffusion MRI can be expressed as the convolution 
of a given response function with the fiber orientation distribution (FOD). The 
FOD is a real-valued function on the unit sphere (S 2 ) giving the orientation and 
the volume fraction of the fiber populations present in a voxel. The response 
function, or kernel, describes the MRI signal attenuation generated by an iso- 
lated single fiber population; it can be estimated from the data or represented 
by means of parametric functions. 

SD approaches represented a big step in reducing the acquisition time of dif- 
fusion MRI, but are known to suffer heavily from noise and intrinsic instabilities 
in solving the deconvolution problem. For this reason, a regularization scheme 
is normally employed. A variety of approaches have been proposed, which are 
generally based on the two assumptions that the FOD is (i) a non-negative func- 
tion and (ii) sparse, i.e. with only a few nonzero values. In fact, at the imaging 
resolution available nowadays, diffusion MRI is sensitive only to the major fiber 
bundles and it is commonly accepted that it can reliably disentangle up to 2-3 
different fiber populations inside a voxel (Jeurissen et al., 2010). Hence, the 
FOD can be considered sparse in nature. 

The recent advent of compressed sensing (CS) theory (Donoho, 2006; Candes 
et al., 2006; Baraniuk, 2007) provided a mathematical framework for the re- 
construction of sparse signal from under-sampled measurements mainly in the 
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context of convex optimization. CS has inspired new advanced approaches in 
the last few years for solving the reconstruction problem in diffusion MRI and 
allowed a further dramatic reduction in the number of samples needed to ac- 
curately infer the fiber structure in each voxel. For instance, Tristan- Vega and 
Westin (2011) and Michailovich et al. (2011) recovered the orientation distribu- 
tion function (ODF) by using different representations for the response function, 
while Mcrlct et al. (2011) and Rathi et al. (2011) focused on the full ensemble 
average propagator (EAP) of the diffusion process. In this work, however, we 
focus on spherical deconvolution based methods and the quantity of interest is 
the FOD. In general, these methods are based on l\ minimization schemes and 
the common goal is to recover the FOD with fewest non-zeros that is compat- 
ible with the acquired MRI data (Ramirez-Manzanares et al., 2007; Pu et al., 
2011; Landman et al., 2012; Mani et al., 2012). However, we believe that the 
l\ norm, identifying the sum of the values of a function, is not adequate for 
characterizing the actual sparsity lying in the FOD, as the sum of the volume 
fractions of the compartments inside a voxel is intrisically equal to 1. Strictly 
speaking, the FOD sparsity is the number of fiber populations, thus identified 
by the lo norm of the FOD. l§ norm problems are generally intractable as they 
are non convex, which explains the usual convex l\ norm relaxation in the con- 
text of compressed sensing. However, among other approaches, a reweighted l\ 
minimization scheme was developed in order to approach £ minimization by a 
sequence of convex weighted i\ problems (Candes et al., 2008). 

In this paper, we propose to exploit the versatility of compressed sensing and 
convex optimization to exploit sparsity more optimally than the state of the art 
and improve the quality of reconstructions in diffusion MRI. To this aim, we 
define a new convex minimization problem through a constrained formulation 
between sparsity prior and data, also making use of a reweighting scheme. We 
evaluate the effectiveness of our proposal for improving the FOD reconstructions 
by comparing with existing state-of-the-art spherical deconvolution methods. 
We report results on both synthetic and real data. 

2. Materials and methods 

2.1. Intra-voxel structure recovery via spherical deconvolution 

As shown by Jian and Vemuri (2007), spherical deconvolution methods can 
be cast into the following computational framework: 



where / is the FOD to be estimated, represents the response function ro- 
tated in direction q £ S 2 and the integration is performed over the unit sphere 
with p = {(f), 9) E S 2 and dil = sm<f>d(f>d8. S(b,q) represents the MRI signal 
measured on the q-space shell acquired with b- value b in direction q 6 S 2 , while 
So is the signal acquired without diffusion weighting. The measurement process 
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can thus be expressed in terms of the general formulation: 

y = <S>x + t), (2) 

where x £ R™ is the FOD, y £ R m is measurement vector, $ £ j^™*" j s ^ e 
observation matrix and rj represents the acquisition noise. 

In the original formulation of Tournier et al. (2004), the FOD x and the 
measurements y were expressed by means of spherical harmonics (SH), and 
the deconvolution problem was solved by a simple matrix inversion. To reduce 
noise artifacts, a low-pass filter was applied for attenuating the high harmonic 
frequencies. The method was improved in (Tournier et al., 2007) by reformulat- 
ing the problem as an iterative procedure where, at each iteration, the current 
solution x^> is used to drive to zero the negative amplitudes of the FOD at the 
next iteration with a Tikhonov regularization (Tikhonov and Arscnin, 1977): 

x {t+1) = argmin \\<5>x - y\\j + A 2 ||i (t) x|||, (3) 

X 

where the free parameter A controls the degree of regularization and can be 
understood as a simple binary mask preserving only the directions of small, prob- 
ably spurious, values of The £2 norm regularization term therefore tends to 
send these values to zero, hence only preserving large values. Interestingly, be- 
yond sending negative values to zero, the operator L thus implicitly promotes a 
weak version of sparsity. || • || p are the usual £ p norms in R™. However, the use of 
the £2 norm to model the sparsity prior does not fully guarantee either positivity 
or sparsity in the recovered FOD. In (Alexander, 2005) a maximum-entropy reg- 
ularization was proposed to recover the FOD as the function that exhibits the 
minimum information content. The method showed higher robustness to noise 
than previous approaches, but was limited by the very high computational cost 
and did not promote sparsity. Other regularization schemes have been proposed 
in the literature, but FOD sparsity has never been addressed with a rigorous 
mathematical formulation. 

Compressed sensing provides a powerful mathematical framework for the 
reconstruction of sparse signals from a low number of data (Donoho, 2006; 
Candcs et al., 2006), mainly in the context of convex optimization. According 
to this theory, it is possible to recover a signal from less samples than the number 
required by the Nyquist sampling theorem, provided that the signal is sparse 
in some sparsity basis \P. Let x £ R™ be the signal to be recovered from the 
m <C n linear measurements y — $x £ R m and a £ R" a sparse representation 
of x through \& £ M. nxn . If the observations y are corrupted by noise and $ 
obeys some randomness and incoherence conditions, then the signal x = 
can be recovered by solving the following convex £\ optimization problem: 

argmin ||a||i subject to ||$Wa — 3/ 1 1 2 < e, (4) 

a 

where e is a bound on the noise level. It is worth noting that, in the context 
of SD reconstructions, the sparsity basis \& boils down to the identity matrix, 
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thus x = a. In (Ramirez-Manzanares et al., 2007; Jian and Vemuri, 2007) the 
sensing basis $, also called dictionary in this context, is generated by applying a 
set of rotations to a given Gaussian kernel (i.e. diffusion tensor) and the sparsest 
coefficients x of this linear combination best matching the measurements y are 
recovered by solving the following constrained minimization problem: 

argmin ||x||i subject to 1 1$ x — y\ | 2 < e, (5) 

cc>0 

where the positivity constraint on the FOD values was directly embedded in the 
formulation of the convex problem. The parameter e can be easily estimated 
from the data assuming that the samples are contaminated with Gaussian noise, 
which is normally the case with MRI data with sufficient signal-to-noise ratio 
(SNR), as shown by Gudbjartsson and Patz (1995). For SNR > 2, in fact, the 
noise in magnitude MRI images can be assume Gaussian-distributed. Conse- 
quently, from a statistical standpoint, the £2 norm is a suitable choice in the 
minimization problems to describe the residual noise on the estimate of x as 
it identifies with the negative log-likelihood. For simplicity one can also ex- 
trapolate the choice of the £2 norm at lower SNR, independently of statistical 
considerations. 

In this context, however, a statistical estimation of e is not reliable. In 
fact, since the aim of these reconstruction methods is to reduce the number of 
measurements, the x 2 estimator of the noise has a large variance around its 
mean. This fact is related to the well-known phenomenon of concentration of 
measure in statistics. To overcome this issue, the reconstruction problem can 
be re-formulated as a regularized (as opposed to constrained) l\ minimization 
as in (Landman ct al, 2012; Pu et al., 2011): 

argmin ||$SB-y||| + \\x\\ u (6) 

where the free parameter f3 controls the trade-off between the data and the 
sparsity constraints. In general, j3 depends on the acquisition scheme and the 
noise level and it must be empirically optimized. Noteworthy, the cost function 
in (5) and (6) aims at minimizing the £\ norm of the FOD x, even though 
J2i x i = \\ x \\i — 1 by construction. For this reason, we reckon that these 
regularized £\ formulations are intrinsically suboptimal and, in this work, our 
main goal is to adequately characterize the actual sparsity lying in the FOD. 

2.2. Reweighted constrained £\ minimization 

In the aim of adequately characterizing the FOD sparsity, we re-formulate 
the spherical deconvolution problem as a constrained £q minimization problem: 

argmin | \<frx — y|| 2 subject to ||as||o < k, (7) 

x>0 

where || • ||o explicitly counts the number of nonzero coefficients and the param- 
eter k represents the expected number of fiber populations in a voxel. 
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As already stated, the £q problems as such are intractable. The reweight- 
ing scheme proposed in Candes et al. (2008) proceeds by sequentially solving 
weighted l\ problems of the form (7), where the £q norm is substituted by a 
weighted l\ norm defined as ||u>a||i = J2i w i \ a i\i f° r positive weights Wi and 
where i indexes vector components. At each iteration, the weights are essen- 
tially set as the inverse of the values of the solution of the previous problem, 
i.e. ie$ ~ l/a?i ■ At convergence, this set of weights makes the weighted £\ 
norm independent of the precise value of the nonzero components, thus mim- 
icking the £ norm while preserving the tractability of the problem with convex 
optimization tools. Of course, it is not possible to have infinite weights where 
the estimated signal values are zero, so that a stability parameter t must also 
be added to the signal value in the selection of the weights. 

The main steps of the reweighted scheme are reported in the algorithm 1. We 
empirically set r = 10~ 5 and the procedure was stopped if ^ ^-i)^ — < 10~ 3 
between two successive iterations or after 20 iterations. At the first iteration 
the weighted l\ norm is the standard £\ norm given w — 1, and therefore the 
constraint ||«/°)a;||i < k is a weak bound on the sum of the fiber compartments 
and does not constitute then a crucial limitation in the procedure. Moreover, 
the algorithm was found to be quite robust to the choice of k, and differences 
were not observed for values up to k — 5. However, as discussed before, we 
can expect to have at maximum 2-3 fiber compartments in each voxel and thus 
in our experiments we fixed k — 3. Finally, an explicit contraint J2i x i = 1 
might have been added. For the sake of simplicity of the algorithm, in this work 
this constaint was not included, assuming it is carried over by well-designed 
bases and the data, as pointed out by Ramirez-Manzanares et al. (2007) . In the 
remaining of the manuscript we will refer to this problem as RSD (Reweighted 
Sparse Deconvolution) . 



Algorithm 1 Reweighted l\ minimization for intra-voxel structure recovery 
Input: Diffusion MRI signal y G K m and sensing basis <I> € U mx ™ 
Output: FOD x e R n 
Set the initial status: 

t and <— 1, i = 1, . . . ,n 

repeat 

Solve the problem: 

argmin \\$>x — subject to ||t«Wa;||i < k 

x>0 

Update the weights: 



t <r- t + 1 
until stopping criterion is satisfied 
x a;(* _1 ) 
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2.3. Comparison framework 

We compared our ^ -based regularizaton scheme against other well estab- 
lished state-of-the-art approaches. In particular, we included in the comparison 
the two regularization schemes (3) and (6), which are based on £2 and l\ pri- 
ors, respectively. In the following, we will refer to these two problems as CSD 
for the former and L2L1 for the latter. To run CSD reconstructions we made 
use of the original mrtrix 1 implementation, setting the optimal parameters as 
suggested by the software itself. To solve the L2L1 and RSD problems we used 
the SPArse Modeling Software (SPAMS) 2 , an open-source and very fast opti- 
mization toolbox written in CH — h for solving various sparse recovery problems. 
Numerical simulations on synthetic data were performed to quantitatively assess 
the performance of CSD, L2L1 and RSD under controlled conditions. A qualitative 
evaluation of the effectiveness on real human brain data was also carried out. 

2.4. Numerical simulations 

Independent voxels with two fiber populations crossing at specific angles 
(30° — 90° range) and with equal volume fractions were synthetically generated. 
The MFJ signal attenuation S corresponding to each voxel configuration was 
simulated by using the exact expression for particles diffusing in a restricted 
cylindrical geometry given in (Soderman and Jonsson, 1995). The following 
parameters were used (Ozarslan ct al., 2006; Jian and Vemuri, 2007): L = 5 mm, 
p = 5 fim, D = 2.02 x 10~ 3 mm 2 /s, A = 20.8 ms, S = 2.4 ms. The signal S 
was contaminated with Rician noise (Gudbjartsson and Patz, 1995) as follows: 

S noisy = V(S + £l) 2 + (6) 2 , (8) 

where £1, £2 ~ A/"(0, a 2 ) and a — 5o/SNR corresponds to a given signal-to-noise 
ratio. We assumed So = 1 without loss of generality. 

For each voxel configuration, the signal was simulated at different b-values, 
b € {500, 1000, 4000} s/mm 2 , and seven q-space sampling schemes were 
tested, respectively with 6, 10, 15, 20, 25, 30 and 50 samples equally distributed 
on half the unit sphere using electrostatic repulsion (Jones et al., 1999) assuming 
antipodal symmetry in diffusion signal. Six different noise levels were consid- 
ered, SNR = 5, 10, ... , 30. For every SNR, 100 repetitions of the same voxel 
were generated using different realizations of the noise. In our experiments, the 
actual signal-to-noise ratio in the simulated signal was always in a range where 
the Gaussian assumption on the noise holds. In the extreme setting with a 
SNR = 5 on the So and b = 4000 s/mm 2 the actual signal-to- noise ratio in the 
signal was about 1.4. 



1 http: / /www. brain. org.au / software / mrtrix 
2 http: / /spams-devel. gforge.inria.fr 



7 



2.5. Evaluation criteria 

As the focus of this work is to improve SD reconstructions, we adopted 
some standard metrics widely used in the literature (Ramircz-Manzanares et al., 
2008; Landman ct al., 2012; Michailovich et al., 2011) to assess the quality of 
reconstructions with respect to number and orientation of the fiber populations 
identified: 

• Probability of false fiber detection. This metric quantifies the correct as- 
sessment of the real number M of populations inside a voxel: 



where M is the estimated number of compartments. As Pd does not 
distinguish between not recovered fibers and extra compartments found 
by the reconstruction, we also make use of the following two quantities 
where needed, n _ and n + , explicitly counting the number under- and 
over-estimated compartments, respectively. 

• Angular error. This metric quantifies the angular accuracy in the estima- 
tion of the directions of the fiber populations in a voxel: 



where d is a true direction and d is its closest estimate. The final value of 
the metric is an average over all the fiber compartments within a voxel. 

Peaks detection was performed using a local maxima search algorithm on the 
FOD recovered with the three methods. For every direction a neighbourhood of 
orientations within a cone of 15° was considered and values smaller than 10% of 
the largest peak were discarded. In the case of CSD, however, we had to increase 
this threshold to 20% as suggested by Tournier et al. (2007) for removing all 
the spurious peaks in order to compare with the other two methods. 

2.6. Real data 

Human brain data was acquired on two young healthy volunteers. Im- 
ages were collected on a 3T Magnetom Trio system (Siemens, Erlangen, Ger- 
many) equipped with a 32-channel head coil using standard DTI protocols 
routinely used in clinical practice. The two datasets have been acquired at 
b = 1000 s/mm 2 using 30 and 20 diffusion gradient directions, respectively, uni- 
formly distributed on half the unit sphere, and in the following they will be 
referred as data 30 and data 2 o- Other acquisition parameters were as follows: 
TR/TE = 7000/82 ms and spatial resolution = 2.5 x 2.5 x 2.5 mm for dataset 
data 30 , while TR/TE = 6000/99 ms and spatial resolution = 2.2 x 2.2 x 3 mm 
for dataset data 2 o- The actual SNR in the 6 = images, computed as the ratio 
of the mean value in a region-of-interest placed in the white matter and the 
standard deviation of the noise estimated in the background, was about 60 in 
data 30 and 30 in data 2 o- 



Pd 



\M ~ M\ 
M 



• 100%, 



(9) 




(10) 
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2.7. Implementation details 

In all our experiments, the response function was estimated from the data 
following the procedure described in Tournier et al. (2004, 2007). In real data, 
the 300 voxels with the highest fractional anisotropy were selected and a dif- 
fusion tensor was fitted from the signal profile in each voxel. In the case of 
numerical simulations, we generated an additional set of data containing voxels 
with a single fiber population and applied the same procedure. A different kernel 
has been estimated for every combination of experimental conditions (b-value, 
SNR, number of samples). 

The same estimated response function was used in all three reconstruction 
methods considered in this work. This was directly used in the CSD algorithm 
whereas, in the case of L2L1 and RSD, the dictionary $ was generated by ro- 
tating it along 200 orientations uniformly distributed on half the unit sphere 
using electrostatic repulsion (Jones et al., 1999). In real data experiments, an 
additional isotropic compartment was included in $ to properly model any par- 
tial/full contamination with CSF. 

Although there are no free parameters to set in CSD and RSD problems (in 
the latter case we can in fact assume k = 3), /3 has to be empirically estimated 
in the case of L2L1 reconstructions. We optimised this parameter following the 
guidelines of Landman et al. (2012), in order to place the method in its best 
conditions. In numerical simulations, we created an additional training dataset 
for every combination of experimental conditions (b-value, SNR, number of 
samples) and 50 reconstructions were performed varying the parameter (3 from 
10~ 4 (i* to /?*, with /3* = ||2$ T y|| 00 computed independently in each voxel. 
The value providing the best reconstruction (using the above metrics) was then 
used to run L2L1 on the actual data used for the final comparison. In real data, 
since the ground-truth is unknown, we set j3 = 0.1 * /3* independently in each 
voxel as suggested in the same work. 

3. Results and discussion 

We quantitatively compared the three reconstruction methods on synthetic 
data under controlled conditions. In particular, the quality of the reconstruc- 
tions was evaluated using the metrics introduced above and selectively varying 
(i) the number of samples and (ii) the b-value of the acquisition scheme, (iii) the 
noise level and (iv) the crossing angle between the fiber compartments. Results 
are reported independently for each experimental condition. 

Figure 1 reports the performances of the three reconstruction methods as 
the number of samples changes. We considered seven acquisition schemes from 
6 to 50 samples and results are reported for a standard acquisition setup, specif- 
ically a shell at b = 2000 s/mm 2 with a SNR = 25. The quality metrics are 
reported here as the average value computed over all simulated crossing angles 
(30° — 90°). Looking at plots the effects of reweighting are clear: RSD always 
outperforms both CSD and L2L1 in identifying the correct number of fiber pop- 
ulations in a voxel (Pd) and results are consistent for all the number of samples 
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G 10 15 20 25 30 50 6 10 15 20 25 30 5C 

number of samples number of samples 




10 15 20 25 30 50 6 10 15 20 25 30 5C 



number of samples number of samples 

Figure 1: Quantitative comparison as a function of the number of samples. The 

values of the four quality metrics are reported for CSD (blue), L2L1 (red) and RSD (green) as the 
number of samples changes. The same color-coding is used consistently in all the figures to 
easily identify each reconstruction method. Values shown here correspond to an experimental 
setting with b = 2000 s/mm 2 and SNR = 25. 



considered. The main benefit of RSD seems to be the drastically decreased num- 
ber of missed fibers (smaller n~), even though also the number of over-estimated 
compartments (n + ) is reduced significantly using RSD as well. Concerning the 
angular accuracy of the recovered fiber directions (eg), reconstructions with RSD 
always resulted in smaller errors as compared to both CSD and L2L1. For 30 or 
more samples the difference is not as significant as it is for Pd- For less than 30 
samples, however, both L2L1 and RSD show a substantial overall improvement 
(about 10° to 15°) with respect to CSD, with RSD further improving the accuracy 
of L2L1 by 1° to 3° on average. 

As expected, the accuracy of the reconstructions degraded as the number of 
samples decreases for all the three methods, although CSD shows a significant 
degradation with less than 30 samples. This can be explained with the SH rep- 
resentation used internally by CSD. In fact, even though the FOD is a function 
on the sphere containing high-resolution features by definition, a maximum SH 
order l max = 4 (or less) can be used for acquisitions with less than 30 sam- 
ples, hence drastically reducing the intrinsic angular resolution of the recovered 
FOD. At least 30 to 60 samples are normally advised for using CSD, so in our 
experiments we have actually tested CSD beyond its applicability range. On the 
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contrary, L2L1 and RSD do not make use of SH and the reconstruction quality 
degrades more smoothly with the under-sampling rate of the MRI signal. In 
the following, then, we will focus on two acquisition schemes to further analyze 
the performances of three methods: (i) in a normal setting with 30 samples and 
(ii) in a condition of high under-sampling with only 15 samples. 

In figure 2 the performances of CSD, L2L1 and RSD are plotted in detail 
as a function of the crossing angle between the fiber populations. Results are 
reported for two acquisition schemes using 30 and 15 samples, both simulated 
at b = 2000 s/mm 2 and SNR = 25. With 30 samples, no modeling errors occur 
(as we pointed out before) and the major source of errors for CSD and L2L1 is 
represented by under-estimation (n~), with both methods starting to severely 
miss fibers for crossing angles below 50°/60°, on the contrary of RSD. In this 
setting, both CSD and L2L1 tend to recover one single peak lying between the two 
real fiber compartments and the maximum angular error for the sole estimated 
compartment is generally upper bounded by half the crossing angle. For this 
reason the overall performances of CSD and L2L1 do not differ significantly from 
RSD despite the drastic improvement in terms of Pd, n~ and n + . Nonetheless, the 
ability of reliably recovering fibers also for crossing angles below 50° — 60° is very 
important for fiber-tracking applications which rely on the accurate estimates 
of white matter orientation, because the propagation of these errors might have 
dramatic consequences in the final estimated trajectories. In an under-sampling 
scenario with 15 samples, however, both CSD and L2L1 are less robust and exhibit 
a higher Pd with a stronger tendency to over-estimate compartments, usually 
in completely arbitrary directions not even close to the true fiber directions. 
The overall improvement of RSD on angular accuracy is more evident, with an 
average enhancement of 1° to 5° with respect to L2L1, while CSD exhibits a 
severe drop of the performances mainly due to modeling limitations. 

Another point worth noting is that CSD achieves slightly better angular ac- 
curacy for crossing angles above 60°. The precision of both L2L1 and RSD is 
in fact affected by the resolution of the grid used to construct the dictionary. 
In our experiments, 200 directions equally distributed on the sphere were used 
to build resulting in a grid resolution of about 10° and thus an intrinsic 
average error of about 5°. To improve the precision it would be sufficient to 
increase the number of directions of the discretization which, however, would 
have serious consequences on the efficiency and stability of the minimization 
algorithm. Interestingly, recent works of Tang et al. (2012) and Candes and 
Fernandez- Granda (2012) explored a novel theory of compressed sensing with 
continuous dictionaries, in the context of which FOD peaks could be thought 
to be located with infinite precision. This topic will be the subject of future 
research. 

So far CSD, L2L1 and RSD have been compared for given acquisition schemes 
and a fixed b = 2000 s/mm 2 . Figure 3 reports the quality of the reconstructions 
with the three approaches as a function of the b-value. The results are shown 
for 30 and 15 samples with a SNR = 25. CSD tends to miss compartments for 
low b-values and over-estimate them at higher 6 (n + and nT are not shown 
here for brevity). This is even more apparent when decreasing the number of 
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crossing angle crossing angle 
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crossing angle crossing angle 



15 samples 
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crossing angle crossing angle 




30 40 50 GO 70 80 90 30 40 50 GO 70 80 90 



crossing angle crossing angle 

Figure 2: Quantitative comparison as a function of the crossing angle. The perfor- 
mances of the three reconstruction methods are detailed separately for each crossing angle used 
in the simulations. Results are reported for 30 and 15 samples, using the same experimental 
configuration of Fig. 1, i.e. b = 2000 s/mm 2 and SNR = 25. 
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Figure 3: Quantitative comparison as a function of the b-value. The dependence of 
the reconstruction quality on the b-values used in the acquisition is reported here for 30 and 
15 samples with a SNR = 25. 



samples in the acquisition to 15, where CSD estimates a lot of spurious peaks at 
high b-values (high n + ) and thus the angular accuracy of the estimated fiber 
directions drops considerably. Interestingly, L2L1 shows the opposite behavior, 
under-estimating at high b and over-estimating at low b, although at a smaller 
rate thus preventing the performances to degrade significantly. Again, in com- 
parison, RSD shows a very stable estimation of the number of fibers. 
Concerning the angular accuracy, all methods showed a minimum for e$ corre- 
sponding tofew 1500 — 2500 s/mm 2 , representing a sort of trade-off between the 
loss in angular resolution happening at small b-values and the stronger noise 
influence at higher b. In fact, as in this work we report the noise level as the 
SNR of the dataset, images at high b-values will have lower actual SNR, and 
thus the noise effects will be inherently stronger. Overall, RSD always results 
in smaller angular errors for the identified fiber compartments than the other 
two methods. The improvement with respect to L2L1 is moderate (1° to 3° 
on average), while the difference with CSD is more emphasized as the b-value 
increases. 

Finally, figure 4 compares the robustness to noise of the three reconstruction 
methods. Six noise levels have been considered, with the SNR of the So dataset 
varying from 5 to 30. The comparison is reported for 30 and 15 samples at 
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b = 2000 s/mm 2 . The results show that RSD clearly outclasses the other two 
methods concerning the estimation of the number of compartments {Pd) and 
results are consistent as the SNR changes. In terms of angular accuracy, eg, 
RSD and L2L1 have very similar eg performances, with RSD performing slightly 
better for SNR > 15 and L2L1 for SNR < 10. CSD, on the contrary, always 
provides significantly lower performance (2° to 6° on average) as compared 
to both RSD and L2L1 at all considered SNRs. The angular accuracy of CSD 
degrades even more in a high under-sampling regime with 15 samples and it 
is almost independent of the noise level. This is again consistent with the 
limitations of the SH representation for acquisitions with very few samples. 

The results of the qualitative evaluation of the three reconstruction methods 
in the case of real data are shown in Fig. 5, for both acquisitions with 30 and 
20 diffusion gradient directions. Subplots A, B and C correspond to the dataset 
data 30 acquired using 30 samples. Even though the acquisition scheme used for 
this dataset is not the setting where our numerical simulations highlighted the 
most substantial differences between the three methods, important conclusions 
can be drawn in favor of RSD. First, the ability of both L2L1 and RSD to properly 
model the isotropic compartment in voxels with full or partial contamination 
with CSF is clearly visible. Moreover, comparing B and C we can observe 
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Figure 5: Qualitative comparison on human brain data. Reconstructions in the corona 
radiata region are shown for: CSD (A and D), L2L1 (B and E) and RSD (C and F). Subplots A, 
B and C correspond to MRI images acquired using 30 samples, while D. E and F are relative 
to the acquisition with 20 samples. All images have been acquired at b = 1000 s/mm 2 . 
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that RSD successfully differentiates gray matter (light gray regions) from CSF 
voxels with pure isotropic and fast diffusion (very bright areas), while L2L1 
appears unable to distinguish them. The yellow frames in the figures highlight 
the corona radiata, a well-known region in the white matter containing crossing 
fibers. As expected from our simulations at this still relatively high number of 
samples, differences are not obvious between the three methods. However, we 
observe that RSD clearly results in sharper and more defined profiles than L2L1, 
whereas the improvements with respect to CSD are confined only to few voxels. 
The not so good performance of L2L1 might be related to the value chosen 
for /?. In contrast, no free parameter has to be empirically optimized in our 
approach. When decreasing the acquisition samples to 20 (subplots D, E and F 
corresponding to data 2 o dataset), reconstructions using RSD are definitely much 
better resolved than both CSD and L2L1. In fact CSD clearly breaks, missing 
many fiber compartments probably due to the aforementioned representation 
limitations. The same happens to L2L1, whose reconstructions appear very 
blurred and noisy. 

4. Conclusion 

We have proposed to exploit the versatility of compressed sensing and con- 
vex optimization to improve spherical deconvolution reconstructions for recov- 
ering the white matter intra-voxel configuration with diffusion MM. We have 
compared our approach with other well-known state-of-the-art regularization 
schemes commonly used in this field on both synthetic and real human brain 
data. Results showed that actually our proposed reweighted constrained l\ reg- 
ularization scheme is very effective and indeed improves the reconstructions. 
This evolution is most remarkable in a high q-space under-sampling regime, 
thus driving the acquisition cost of HARDI always closer to that of DTI. 

Even though we focused here on single voxel experiments, future work will 
be devoted to study the applicability and the effectiveness of our approach in 
more sophisticated frameworks exploiting the spatial correlation in the data. 
Additionally, future research will investigate the use of the recently proposed 
continuous CS theory with the aim of further improving the accuracy of the 
reconstructions and reducing the acquisition time. 
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